// MtxDiff
// Loads two matrices in matrix market format.
// For every non-zero element, it computes the lexicographical distance
// between the entries in the two matrices.  Reports if the matrices
// are exactly (or approximately) equal.
//
#include <map>
#include <cmath>
#include <string>
#include <cstdlib>
#include <complex>
#include <fstream>
#include <iomanip>
#include <assert.h>
#include <iostream>
//#include <boost/progress.hpp>
//#include <boost/tokenizer.hpp>
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#include <Types.h>
#include <Defines.h>
#include <LoadMatrix.h>
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
using namespace std;
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
template<typename T> void PrintMappedMtxToScreen( map<Coordinate, T> & A, const string & mtxName ) {
  typename map< Coordinate, T>::iterator it;
  cout << "Printing matrix " << mtxName << endl;
  for( it = A.begin(); it != A.end(); it++ ) {
    Coordinate c = it->first;
    T          v = it->second;
    cout << mtxName << "[" << c.row << ", " << c.col << "] = " << v << endl;
  }
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
double myAbs( complex<double> v ) {
   double rs = v.real() * v.real();
   double is = v.imag() * v.imag();
   return sqrt( rs + is );
}
double myAbs( double v ) {
   return fabs( v );
}
double myAbs( S64 v ) {
   return abs( v );
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void PrintParameters( const U64 maxDist, const double zeroThresh, const bool allowMissingEntries ) {
	string ame = allowMissingEntries ? "true" : "false";
	cout << "Parameters:" << endl;
	cout << "maximum lexicographical distance: " << maxDist << endl;
	cout << "round to zero if |val| < " << zeroThresh << endl;
	cout << "missing entries allowed: " << ame        << endl;	
}
void PrintFailMessage( const U64 maxDist, const double zeroThresh, const bool allowMissingEntries ) {
	PrintParameters( maxDist, zeroThresh, allowMissingEntries );
  cout << "failed..." << endl;
}
void PrintPassMessage( const U64 maxDist, const double zeroThresh, const bool allowMissingEntries ) {
	PrintParameters( maxDist, zeroThresh, allowMissingEntries );
  cout << "passed..." << endl;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
bool CompareMatrixProperties( const MatrixProperties & mtxInfo1, const MatrixProperties & mtxInfo2 ) {
   // returns true if matrix dimensions are equal
   // note: this does not compare nnz
   if( mtxInfo1.m != mtxInfo2.m || mtxInfo1.n != mtxInfo2.n ) {
      cout << endl;
      cout << "The matrices do not match."            << endl;
      cout << "details for " << mtxInfo1.fullFileName << endl;
      cout << "  m = "       << mtxInfo1.m            << endl;
      cout << "  n = "       << mtxInfo1.n            << endl;
      cout << "nnz = "       << mtxInfo1.nnz          << endl;
      cout << "details for " << mtxInfo2.fullFileName << endl;
      cout << "  m = "       << mtxInfo2.m            << endl;
      cout << "  n = "       << mtxInfo2.n            << endl;
      cout << "nnz = "       << mtxInfo2.nnz          << endl;
      return false;
   }
   return true;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
template <typename F, typename I> void PrintLexicographicalMidPoint( void ) {
   union valueHolder {
      F f;
      I i;
   };
     
   valueHolder foo;
   foo.f = -0.0;
   cout << "midpt = 0x" << hex << foo.i << dec << endl;
}
bool AlmostEqualRelativeOrAbsolute( double val1, double val2, double maxRelativeError, double maxAbsoluteError, const bool verbose ) {

   // returns true if val1 ~= val2
   // uses a relative error metric to check for near equality
   //
   // credit here to:
   // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

   double A = val1;
   double B = val2;
   
   if( fabs(A - B) < maxAbsoluteError ) {
      return true;
   }
   
   double relativeError;
   
   if( fabs(B) > fabs(A) ) {
      relativeError = fabs( (A - B) / B );
   }
   else {
      relativeError = fabs( (A - B) / A );
   }
   
   if( relativeError <= maxRelativeError ) {
      return true;
   }
   return false;
}
bool LexicographicalCheck( float  val1, float  val2, S32 maxDist, const float  zeroThresh, const bool verbose ) {

   // returns true if val1 ~= val2                                                                                                         
   // uses lexicographical distance between the IEEE encoded floating point values                             
   // parameterized by "midpt" -- the lexicographical midpt for either float or double                         
   //                                                                                                                                      
   // credit here to:                                                                                                                      
   // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm                                

   // Make sure maxDist is non-negative and small enough that the                                              
   // default NAN won't compare as equal to anything.
   DEBUG_PRINT( "in float LexicographicalCheck()" );
   assert( maxDist >= 0 && maxDist < 4 * 1024 * 1024 );

   S32 intDiff;                                                                                    
   const S32 zero = 0LL;

   U32 numBits  = sizeof(S32) * 8;                                                                                                     
   U32 numShift = numBits == 64 ? 28 : 0;                                                                                                  
   maxDist <<= numShift;                                                                                                                   

   union ValueHolder {
      float  f;
      S32    i;
   };

   ValueHolder v1;
   ValueHolder v2;

   v1.f = val1;
   v2.f = val2;

   if( fabs(v1.f) < zeroThresh ) v1.f = 0.0;                                                                                               
   if( fabs(v2.f) < zeroThresh ) v2.f = 0.0;                                                                                               

   // Make ints lexicographically ordered as a twos-complement int
   if( v1.i < zero ) v1.i = 0x80000000 - v1.i;
   if( v2.i < zero ) v2.i = 0x80000000 - v2.i;
   intDiff = abs( v1.i - v2.i );

   bool match = false;                                                                                                                     
   if( intDiff <= maxDist ) {
      match = true;                                                                                                                         
   }
   if( false && verbose && !match ) {
      cout << "      int1 = " << dec << setw(20) <<    v1.i << " = 0x" << setfill('0') << hex << setw(8)  << v1.i << " = " << v1.f << setfill(' ') << endl;
      cout << "      int2 = " << dec << setw(20) <<    v2.i << " = 0x" << setfill('0') << hex << setw(8)  << v2.i << " = " << v2.f << setfill(' ') << endl;
      cout << "      dist = " << dec << setw(20) << intDiff << endl;
   }
   return match;
}                                                                                                                                         
bool LexicographicalCheck( S64    val1, S64    val2, S64 maxDist, const double zeroThresh, const bool verbose ) {
   DEBUG_PRINT( "in S64 LexicographicalCheck()" );
   return val1 == val2;                                                           
}                                                                                
bool LexicographicalCheck( double val1, double val2, S64 maxDist, const double zeroThresh, const bool verbose ) {

   // returns true if val1 ~= val2
   // uses lexicographical distance between the IEEE encoded floating point values
   // parameterized by "midpt" -- the lexicographical midpt for either float or double
   //
   // credit here to:
   // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

   // Make sure maxDist is non-negative and small enough that the
   // default NAN won't compare as equal to anything.
   DEBUG_PRINT( "in double LexicographicalCheck()" );
   assert( maxDist >= 0 && maxDist < 4 * 1024 * 1024 );

   S64 intDiff;
   const S64 zero = 0LL;

   U32 numBits  = sizeof(S64) * 8;
   U32 numShift = numBits == 64 ? 28 : 0;
   maxDist <<= numShift;
  
   union ValueHolder {
      double f;
      S64    i;
   };

   ValueHolder v1;
   ValueHolder v2;
   
   v1.f = val1;
   v2.f = val2;
   
   if( fabs(v1.f) < zeroThresh ) v1.f = 0.0;
   if( fabs(v2.f) < zeroThresh ) v2.f = 0.0;
   
   // Make ints lexicographically ordered as a twos-complement int
   if( v1.i < zero ) v1.i = 0x8000000000000000LL - v1.i;
   if( v2.i < zero ) v2.i = 0x8000000000000000LL - v2.i;
   intDiff = labs( v1.i - v2.i );
   
   bool match = false;
   if( intDiff <= maxDist ) {
      match = true;
   }
   if( false && verbose && !match ) {
      S64 adjDiff = intDiff >> numShift;
      cout << "      int1 = " << dec << setw(20) <<    v1.i << " = 0x" << setfill('0') << hex << setw(16) << v1.i << " = " << v1.f << setfill(' ') << endl;
      cout << "      int2 = " << dec << setw(20) <<    v2.i << " = 0x" << setfill('0') << hex << setw(16) << v2.i << " = " << v2.f << setfill(' ') << endl;
      cout << "      dist = " << dec << setw(20) << intDiff << endl;
      cout << "   adjDist = " << dec << setw(20) << adjDiff << endl;
   }
   return match;
}
bool LexicographicalCheck( complex<double> val1, complex<double> val2, S64 maxDist, const double zeroThresh, const bool verbose) {
   DEBUG_PRINT( "in complex LexicographicalCheck()" );
   const double r1 = val1.real();
   const double r2 = val2.real();
   const double c1 = val1.imag();
   const double c2 = val2.imag();
   bool matched = true;
   if( !LexicographicalCheck( r1, r2, maxDist, zeroThresh, verbose ) ) {
   	cout << "... comparing real parts." << endl;
   	matched = false;
   }
   if( !LexicographicalCheck( c1, c2, maxDist, zeroThresh, verbose ) ) {
   	cout << "... comparing imag parts." << endl;
   	matched = false;
   }
   return matched;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
template <typename T> void PrintVector(const vector<T> & v, const string & name) {
   for( U64 i=0; i < v.size(); i++ ) {
      cout << name << "[" << i << "] = " << v[i] << endl;
   }  
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
template<typename T> double dround( T f, int prec ) {
   char buf[100];
   char fmt[20];
   sprintf(fmt, "%%.%de", prec);
   sprintf(buf, fmt, f);
   double v = atof( buf );
   return v;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
template<typename T, typename F, typename I>
bool MappedMatrixDiff( map<Coordinate, T> & A, map<Coordinate, T> & B, const MatrixProperties & mtxInfoA, const MatrixProperties & mtxInfoB, const I & maxDist, const F & zeroThresh, const bool allowMissingEntries, const bool verbose ) {

   bool nearEqual = true;
	
   Coordinate c;
   map< Coordinate, bool > pattern;
   map< Coordinate, bool >::iterator it;
   typename map< Coordinate, T >::iterator itrB;
   typename map< Coordinate, T >::iterator itrA;

   unsigned int nchars = MAX( mtxInfoA.fullFileName.size(), mtxInfoB.fullFileName.size() );

   if( !CompareMatrixProperties( mtxInfoA, mtxInfoB ) ) {
      return false;
   }

   if( verbose ) {
      PrintLexicographicalMidPoint<F, I>();
   }

   DEBUG_PRINT( "Creating a pattern matrix with non-zero entries from both A & B." );
   for( itrA = A.begin(); itrA != A.end(); itrA++ ) {
      c.row = itrA->first.row;
      c.col = itrA->first.col;
      pattern[c] = true;
   }
   for( itrB = B.begin(); itrB != B.end(); itrB++ ) {
      c.row = itrB->first.row;
      c.col = itrB->first.col;
      pattern[c] = true;
   }

   DEBUG_PRINT( "Iterating over the non-zero entries in the pattern matrix." );
   for( it = pattern.begin(); it != pattern.end(); it++ ) {

      c.row = it->first.row;
      c.col = it->first.col;
      DEBUG_PRINT( "Found entry in pattern at [" << c.row << ", " << c.col << "]" );
      
      itrA = A.find( c );
      itrB = B.find( c );
      
      bool entryInA = itrA != A.end();
      bool entryInB = itrB != B.end();

      if( entryInA && entryInB ) {
         DEBUG_PRINT( "Entry exists in both A & B." );
         T v1 = itrA->second;
         T v2 = itrB->second;
      
         if( !LexicographicalCheck( v1, v2, maxDist, zeroThresh, verbose) ) {
            // could use: if( ! AlmostEqualRelativeOrAbsolute(val1, val2, 0.05,  0.005, verbose) ) {
            if( verbose ) {
               cout << setw(nchars) << mtxInfoA.fullFileName << "[" << setw(5) << c.row+1 << ", " << setw(5) << c.col+1 << "] = " << v1 << endl;
               cout << setw(nchars) << mtxInfoB.fullFileName << "[" << setw(5) << c.row+1 << ", " << setw(5) << c.col+1 << "] = " << v2 << endl;
            }
            nearEqual = false;
         }
         //double delta = myAbs( v1 - v2 );
         //double pcnt  = 100 * delta / min( myAbs(v1), myAbs(v2) );
         //if( pcnt > thresh ) {
         //  value += delta;
         //}
         ////if( (v1-v2) != 0.0 ) {
         ////  cout << "diff at [" << c.row << ", " << c.col << "] = ";
         ////  cout << fabs(v1-v2) << ", v1 = " << v1 << ", v2 = " << v2 << ", pcnt error = " << 100*fabs(v1-v2) / max(v1, v2) << endl;
         ////}
      }
      else if( entryInA ) {
         DEBUG_PRINT( "Entry exists in A, but not B." );
         T v1 = itrA->second;
         T v2 = v1 - v1;
         if( !allowMissingEntries || !LexicographicalCheck( v1, v2, maxDist, zeroThresh, verbose ) ) {
            //if( ! AlmostEqualRelativeOrAbsolute(0.0, val2, 0.05,  0.005, verbose) ) {
            if( verbose ) {
               cout << setw(nchars) << mtxInfoA.fullFileName << "[" << setw(5) << c.row+1 << ", " << setw(5) << c.col+1 << "] = " << v1   << endl;
               cout << setw(nchars) << mtxInfoB.fullFileName << "[" << setw(5) << c.row+1 << ", " << setw(5) << c.col+1 << "] = no entry" << endl;
            }
            nearEqual = false;
         }
          //value += myAbs( itrA->second );
      }
      else if( entryInB ) {
         DEBUG_PRINT( "Entry exists in B, but not A." );
         T v2 = itrB->second;
         T v1 = v2 - v2;
         if( !allowMissingEntries || !LexicographicalCheck( v1, v2, maxDist, zeroThresh, verbose ) ) {
            //if( ! AlmostEqualRelativeOrAbsolute(0.0, val2, 0.05,  0.005, verbose) ) {
            if( verbose ) {
               cout << setw(nchars) << mtxInfoA.fullFileName << "[" << setw(5) << c.row+1 << ", " << setw(5) << c.col+1 << "] = no entry" << endl;
               cout << setw(nchars) << mtxInfoB.fullFileName << "[" << setw(5) << c.row+1 << ", " << setw(5) << c.col+1 << "] = " << v2   << endl;
            }
            nearEqual = false;
         }
         //value += myAbs( itrA->second );	
      }
      else {
         if( verbose ) {
            cout << "no entry, either mtx at row = " << c.row << " col = " << c.col << endl;
         }
         return false;
      }
  }

  return nearEqual;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
int main( int argc, char** argv ) {

   const bool verbose = true;
   const bool expandSymmMats = false;
   const bool allowMissingEntries = false;
   double zeroThresh = 0.0;
   S64    maxDistIn  = 0LL;
   U64    maxDist    = maxDistIn;
   bool matched = false;

   MatrixProperties mtxInfo1;
   MatrixProperties mtxInfo2;
   
   map< Coordinate, double >          Ad;
   map< Coordinate, double >          Bd;
   map< Coordinate, S64 >             Ai;
   map< Coordinate, S64 >             Bi;
   map< Coordinate, complex<double> > Ac;
   map< Coordinate, complex<double> > Bc;
   
   //----- Check input arguments
   if( argc < 3 ) {
      cout << "Usage: "<< argv[0] << " <matrix1 filename> <matrix2 filename> [maxDist] [zeroThresh]" << endl;
      exit(1);
   }
   
   const string fileName1( argv[1] );
   const string fileName2( argv[2] );
   ifstream inputFile1( fileName1.c_str() );
   ifstream inputFile2( fileName2.c_str() );
   
   if( !inputFile1.is_open() ) {
      cout << "Cannot open matrix file: "<< argv[1] << endl;
      exit( -1 );
   }
   if( !inputFile2.is_open() ) {
      cout << "Cannot open matrix file: "<< argv[2] << endl;
      exit( -1 );
   }
   
   if( argc > 3 ) {
      maxDistIn = atol( argv[3] );
      assert( maxDistIn >= 0LL );
      maxDist = maxDistIn;
   }
   if( argc > 4 ) {
      zeroThresh = atof( argv[4] );
      assert( zeroThresh > 0.0 );
   }
    
   ProcessMMFileHeader( fileName1, inputFile1, mtxInfo1 );
   ProcessMMFileHeader( fileName2, inputFile2, mtxInfo2 );
   
   if( mtxInfo1.element_type == "real" ) {
      cout << "Reading in A from file: " << argv[1] << endl;  ReadInMMElements( Ad, inputFile1, mtxInfo1, expandSymmMats );
      cout << "Reading in B from file: " << argv[2] << endl;  ReadInMMElements( Bd, inputFile2, mtxInfo2, expandSymmMats );
      //PrintMappedMtxToScreen<double>(Ad, "A");
      //PrintMappedMtxToScreen<double>(Bd, "B");
      
      matched = MappedMatrixDiff<double, double, S64>( Ad, Bd, mtxInfo1, mtxInfo2, maxDist, zeroThresh, allowMissingEntries, verbose );
   }
   if( mtxInfo1.element_type == "integer" ) {
      cout << "Reading in A from file: " << argv[1] << endl;
      ReadInMMElements(Ai, inputFile1, mtxInfo1, expandSymmMats);
      cout << "Reading in B from file: " << argv[2] << endl;
      ReadInMMElements(Bi, inputFile2, mtxInfo2, expandSymmMats);
      //PrintMappedMtxToScreen<S64>(Ai, "A");
      //PrintMappedMtxToScreen<S64>(Bi, "B");
      matched = MappedMatrixDiff< S64, double, S64>( Ai, Bi, mtxInfo1, mtxInfo2, maxDist, (S64)zeroThresh, allowMissingEntries, verbose );
   }
   if( mtxInfo1.element_type == "complex" ) {
      cout << "Reading in A from file: " << argv[1] << endl; ReadInMMElements(Ac, inputFile1, mtxInfo1, expandSymmMats);
      cout << "Reading in B from file: " << argv[2] << endl; ReadInMMElements(Bc, inputFile2, mtxInfo2, expandSymmMats);
      //PrintMappedMtxToScreen<complex<double> >(Ac, "A");
      //PrintMappedMtxToScreen<complex<double> >(Bc, "B");
      matched = MappedMatrixDiff< complex<double>, double, S64>( Ac, Bc, mtxInfo1, mtxInfo2, maxDist, zeroThresh, allowMissingEntries, verbose );
   }

   inputFile1.close();
   inputFile2.close();
   
   if( !matched ) {
      PrintFailMessage( maxDist, zeroThresh, allowMissingEntries );
      return -1;
   }
   PrintPassMessage( maxDist, zeroThresh, allowMissingEntries );
   return 0;
}
